#include <cstdint>
#include <cstring>
#include <iostream>
#include <cstdint>
#include <chrono>
#include <random>
#include <algorithm>


void mean_filter_optimized(const uint8_t* input, uint8_t* output, int width, int height) {
    const int size = width * height;
    uint16_t* temp = new uint16_t[size]; // 临时数组存储水平和

    // 水平方向处理
    for (int i = 0; i < height; ++i) {
        const uint8_t* row = input + i * width;
        uint16_t* tmp_row = temp + i * width;

        // 处理第一个像素：左边界填充
        tmp_row[0] = row[0] + row[0] + row[1];

        // 处理中间像素
        for (int j = 1; j < width - 1; ++j) {
            tmp_row[j] = row[j - 1] + row[j] + row[j + 1];
        }

        // 处理最后一个像素：右边界填充
        tmp_row[width - 1] = row[width - 2] + row[width - 1] + row[width - 1];
    }

    // 垂直方向处理
    for (int j = 0; j < width; ++j) {
        // 处理第一个像素：上边界填充
        uint16_t sum = temp[j] + temp[j] + temp[j + width];
        output[j] = static_cast<uint8_t>(sum / 9);

        // 处理中间像素
        for (int i = 1; i < height - 1; ++i) {
            sum = temp[j + (i - 1) * width] + temp[j + i * width] + temp[j + (i + 1) * width];
            output[i * width + j] = static_cast<uint8_t>(sum / 9);
        }

        // 处理最后一个像素：下边界填充
        sum = temp[j + (height - 2) * width] + temp[j + (height - 1) * width] + temp[j + (height - 1) * width];
        output[(height - 1) * width + j] = static_cast<uint8_t>(sum / 9);
    }

    delete[] temp;
}

void integral_image(const uint8_t* input, uint32_t* int_img, int width, int height) {
    // 计算积分图，扩展1行1列方便边界处理
    for (int i = 0; i <= height; ++i) {
        for (int j = 0; j <= width; ++j) {
            if (i == 0 || j == 0) {
                int_img[i * (width + 1) + j] = 0;
            } else {
                int_img[i * (width + 1) + j] = input[(i-1) * width + (j-1)]
                    + int_img[(i-1) * (width + 1) + j]
                    + int_img[i * (width + 1) + (j-1)]
                    - int_img[(i-1) * (width + 1) + (j-1)];
            }
        }
    }
}

void mean_filter_integral(const uint8_t* input, uint8_t* output, int width, int height) {
    const int int_w = width + 1, int_h = height + 1;
    std::vector<uint32_t> int_img(int_w * int_h, 0);
    integral_image(input, int_img.data(), width, height);

    for (int i = 0; i < height; ++i) {
        for (int j = 0; j < width; ++j) {
            // 计算3x3区域在积分图中的坐标
            int x1 = std::max(j - 1, 0), y1 = std::max(i - 1, 0);
            int x2 = std::min(j + 1, width - 1), y2 = std::min(i + 1, height - 1);

            // 积分图区域和公式
            uint32_t sum = int_img[(y2 + 1) * int_w + (x2 + 1)]
                         - int_img[(y2 + 1) * int_w + x1]
                         - int_img[y1 * int_w + (x2 + 1)]
                         + int_img[y1 * int_w + x1];
            output[i * width + j] = static_cast<uint8_t>(sum / 9);
        }
    }
}

// 生成随机测试图像
void generate_random_image(uint8_t* image, int size) {
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_int_distribution<> dis(0, 255);
    for (int i = 0; i < size; ++i) {
        image[i] = static_cast<uint8_t>(dis(gen));
    }
}

// 验证结果一致性（忽略边界1像素）
bool verify_results(const uint8_t* res1, const uint8_t* res2, int width, int height) {
    const int margin = 1; // 忽略边界
    for (int i = margin; i < height - margin; ++i) {
        for (int j = margin; j < width - margin; ++j) {
            if (std::abs(res1[i*width + j] - res2[i*width + j]) > 1) { // 允许1的误差
                std::cerr << "Mismatch at (" << i << "," << j << "): "
                          << static_cast<int>(res1[i*width + j]) << " vs "
                          << static_cast<int>(res2[i*width + j]) << std::endl;
                return false;
            }
        }
    }
    return true;
}

int main() {
    const int W = 1024, H = 1024;
    const int NUM_TRIALS = 10; // 测试次数
    
    // 分配内存
    uint8_t* image = new uint8_t[W*H];
    uint8_t* result_optimized = new uint8_t[W*H]();
    uint8_t* result_integral = new uint8_t[W*H]();

    // 生成测试数据
    generate_random_image(image, W*H);

    // 预热运行（避免冷缓存影响）
    mean_filter_optimized(image, result_optimized, W, H);
    mean_filter_integral(image, result_integral, W, H);

    // 测试行列分离算法
    auto t1 = std::chrono::high_resolution_clock::now();
    for (int i = 0; i < NUM_TRIALS; ++i) {
        mean_filter_optimized(image, result_optimized, W, H);
    }
    auto t2 = std::chrono::high_resolution_clock::now();
    auto optimized_time = std::chrono::duration_cast<std::chrono::milliseconds>(t2 - t1).count();

    // 测试积分图算法
    t1 = std::chrono::high_resolution_clock::now();
    for (int i = 0; i < NUM_TRIALS; ++i) {
        mean_filter_integral(image, result_integral, W, H);
    }
    t2 = std::chrono::high_resolution_clock::now();
    auto integral_time = std::chrono::duration_cast<std::chrono::milliseconds>(t2 - t1).count();

    // 验证结果
    bool match = verify_results(result_optimized, result_integral, W, H);

    // 输出结果
    std::cout << "=== Performance Report (1Kx1K) ===" << "\n"
              << "Optimized (Separable) Time: " << optimized_time << " ms (" 
              << optimized_time/NUM_TRIALS << " ms per run)\n"
              << "Integral Image Time:        " << integral_time << " ms ("
              << integral_time/NUM_TRIALS << " ms per run)\n"
              << "Speedup Ratio:              " << static_cast<float>(integral_time)/optimized_time << "x\n"
              << "Results Match:              " << (match ? "Yes" : "No") 
              << " (边界差异已忽略)" << std::endl;

    // 清理内存
    delete[] image;
    delete[] result_optimized;
    delete[] result_integral;
    return 0;
}

// 示例用法
// int main() {
//     const int W = 1024, H = 1024;
//     uint8_t* image = new uint8_t[W * H];
//     uint8_t* filtered = new uint8_t[W * H];
// 
//     // 假设image已填充数据
//     mean_filter_optimized(image, filtered, W, H);
// 
//     delete[] image;
//     delete[] filtered;
//     return 0;
// }
